## test visualizazion of GreenLand masks


link <- "/data/Dagobah/greengrass/grassland_ger/01_CTM_GL_mask_17-19/GL_mask_2017-2019.tif"

mask_GL <- terra::rast(link)

plot(mask_GL)

TSI I/C/O

Import TSI images (single location)

# Full mosaic
#link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.vrt.ovr"

#  base case
#link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/mosaic/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.vrt.ovr"

# Single tile base case
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/base_case/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.tif"

tsi_base_raw <- rast(link)

# Single tile case day_16_sigma_8_16_32
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

tsi_day_16_sigma_8_16_32_raw <- rast(link)


# Single tile case day_16_sigma_8_8_8_16_16_16_64
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_8_8_16_16_16_64/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

tsi_day_16_sigma_8_8_8_16_16_16_64_raw <- rast(link)

# Single tile case day_16_sigma_8_8_8_16_16_16_32_64
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_8_8_16_16_16_32_64/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

tsi_day_16_sigma_8_8_8_16_16_16_32_64_raw <- rast(link)


# combine TSIs into a list to facilitate functional loop wise data operations

tsi_list <- list(tsi_base_raw, tsi_day_16_sigma_8_16_32_raw, tsi_day_16_sigma_8_8_8_16_16_16_64_raw, tsi_day_16_sigma_8_8_8_16_16_16_32_64_raw)

cases <- c("base_case", "day_16_sigma_8_16_32", "day_16_sigma_8_8_8_16_16_16_64", "day_16_sigma_8_8_8_16_16_16_32_64")
names(tsi_list) <- cases

tsi_obs_length <- list()
lapply(seq_along(tsi_list),
       function(i) 
         tsi_obs_length[i] <<-  length(names(tsi_list[[i]]))
)

test vis of TSI

create stack for direct location comparison

tsi_stack <- c(tsi_list[[1]][[c(430)]],
               tsi_list[[2]][[c(675)]],
               tsi_list[[3]][[c(675)]],
               tsi_list[[4]][[c(675)]])

static maps

#de <- st_read("/data/Dagobah/dc/deu/vector/DEU.gpkg")

#plot(tsi_list[[2]][[800:808]])

lapply(seq_along(tsi_list),
       function(i) plot(tsi_stack[[i]],
                        legend=T, col = viridis(n=100,option="D"),
                        range=c(0,10000), 
                        plg=list(title= paste(names(tsi_list[i]), "\n", names(tsi_stack[[i]]))
                                 ))
)

Interactive map

pal <- colorNumeric(c("#0C2C84", "#41B6C4", "#FFFFCC"), values(tsi_list[[4]][[675]]),
  na.color = "transparent")

leaflet() |>
  addTiles() |>
  addRasterImage(tsi_list[[4]][[675]], group = "show_layer") |> 
  addLegend(pal = pal, 
            values = values(tsi_list[[4]][[600:600]]),
            title = "Reflectance") |> 
  addLayersControl(
    overlayGroups = c("show_layer"),
    options = layersControlOptions(collapsed = FALSE)
    )

data wrangeling of TSI images

Compute and set min and max values for tsi cases

lapply(seq_along(tsi_list),
       function(i) 
      tsi_list[[i]] <<-  setMinMax(tsi_list[[i]])
)
## [[1]]
## class       : SpatRaster 
## dimensions  : 1000, 1000, 675  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4166026, 4196026, 3344920, 3374920  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSS.tif 
## names       : 19840~LND05, 19840~LND05, 19840~LND05, 19840~LND05, 19841~LND05, 19841~LND05, ... 
## min values  :        -507,        1140,        2569,        1768,         716,        -355, ... 
## max values  :        8815,        8873,        8598,        8288,        7954,        9881, ... 
## 
## [[2]]
## class       : SpatRaster 
## dimensions  : 1000, 1000, 866  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4166026, 4196026, 3344920, 3374920  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif 
## names       : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ... 
## min values  :      NaN,      NaN,      NaN,     -507,      475,      225, ... 
## max values  :      NaN,      NaN,      NaN,     8815,     8873,     8873, ... 
## 
## [[3]]
## class       : SpatRaster 
## dimensions  : 1000, 1000, 866  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4166026, 4196026, 3344920, 3374920  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif 
## names       : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ... 
## min values  :      570,      592,      613,      871,      875,      349, ... 
## max values  :     8873,     8873,     8873,     8873,     8873,     8873, ... 
## 
## [[4]]
## class       : SpatRaster 
## dimensions  : 1000, 1000, 866  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : 4166026, 4196026, 3344920, 3374920  (xmin, xmax, ymin, ymax)
## coord. ref. : ETRS89-extended / LAEA Europe (EPSG:3035) 
## source      : 1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif 
## names       : 19840101, 19840117, 19840202, 19840218, 19840306, 19840322, ... 
## min values  :      570,      592,      613,      703,      857,      448, ... 
## max values  :     8873,     8873,     8873,     8873,     8873,     8873, ...
### alternatice slower na layer computation method
# 
# na_index   <- list("base_case"=c(), "day_16_sigma_8_16_32"=c(), "day_16_sigma_4_8_16"=c(), "day_16_sigma_8_8_8_16_16_16_64"=c(), "day_8_sigma_8_16_32"=c())
# 
# system.time(
#   
# lapply(seq_along(tsi_list),
#        function(i) 
#          
#          for (j in c(1:nlyr(tsi_list[[i]]))) {
#            
#            na_index[[i]][[j]] <-  minmax(is.na(foo[[i]][[j]]))[1]
#            
#          }
# )
# 
# )

For all cases most cases, there are occasionally layers consisting exclusively of NAs (ie where min and max is NA). Overall there are relatively rare (~3-15%; computed below), and occur mostly at the beginning of the Landsat program (1984-1986), and occasionally during winter months.

The empty layers are included in the force output to satisfy the specified day window intervals for the interpolation. Unfortunately this has the downstream effect that when attempting to convert the spatial object to a df, R throws errors when empty layers are included. The current workaround is to compute the the min/max per layer, create an index of valid layers, and then use this index to subset the tsi stacks.

remove all NA layers

# create filter out all layers that dont contain any values
tsi_na_index <- list()
lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_na_index[i] <<- minmax(tsi_list[[i]]) |> 
         as.data.frame() |> 
         slice(1) |> 
         pivot_longer(everything()) |> 
         mutate(id = row_number()) |> 
         na.omit()
)

# convert into index vector
lapply(seq_along(tsi_list),
       function(i) 
         tsi_na_index[[i]] <<- as.vector(tsi_na_index[[i]])
       
)

# subset tsi object with only layers containing values
lapply(seq_along(tsi_list),
       function(i) 
         tsi_list[[i]] <<- tsi_list[[i]][[ tsi_na_index[[i]] ]] 
)


# filter out reflectance values below and above 0-100 % range
# lapply(seq_along(tsi_list),
#        function(i) 
#          tsi_list[[i]] <- tsi_list[[i]] |> 
#   filter(across(1:length(names(tsi_list[[i]])), ~ . > 0),
#          across(1:length(names(tsi_list[[i]])), ~ . < 10000))
# )

Check NA layer occurance

# lapply(seq_along(tsi_list),
#        function(i) 
#        tsi_obs_length[[i]] - length(names(tsi_list[[i]]))
# )

lapply(seq_along(tsi_list),
       function(i) 
         (paste("In the", names(tsi_list)[[i]],"scenario", (1-round(( length(names(tsi_list[[i]]))/tsi_obs_length[[i]]), digits=4))*100, "% of layers consisted completely of NA values (n=", ( tsi_obs_length[[i]] - length(names(tsi_list[[i]]))), "out of",  tsi_obs_length[[i]], ")")) 
       
)

[[1]] [1] “In the base_case scenario 0.149999999999995 % of layers consisted completely of NA values (n= 1 out of 675 )”

[[2]] [1] “In the day_16_sigma_8_16_32 scenario 3.35 % of layers consisted completely of NA values (n= 29 out of 866 )”

[[3]] [1] “In the day_16_sigma_8_8_8_16_16_16_64 scenario 0 % of layers consisted completely of NA values (n= 0 out of 866 )”

[[4]] [1] “In the day_16_sigma_8_8_8_16_16_16_32_64 scenario 0 % of layers consisted completely of NA values (n= 0 out of 866 )”

create updated stack for direct location comparison

(ie because NA filter removed some layers and dates dont correspond when using the previous index)

tsi_stack <- c(tsi_list[[1]][[c(430)]],
               tsi_list[[2]][[c(648)]],
               tsi_list[[3]][[c(675)]],
               tsi_list[[4]][[c(675)]])
#rasterVis::levelplot(tsi[[c(24)]], par.settings = viridisTheme)
#rasterVis::levelplot(tsi[[c(1,4,24,25:35)]], par.settings = viridisTheme)


lapply(seq_along(tsi_list),
       function(i) 
         plot(tsi_stack[[i]], 
              legend=T, col = viridis(n=100,option="D"),
              range=c(0,10000), 
              plg=list(title= paste(names(tsi_list[i]), "\n", names(tsi_list[[i]][[674]]) ))
         )
)

generation of sample points (performed in Rstudio GUI)

roi_names <- c("watt_GL", "large", "medium", "small")

# select the ROI for each location type
# watt_GL   <- terra::sel(tsi_list[[4]][[c(430:430)]])
# large  <- terra::sel(tsi_list[[4]][[c(430:430)]])
# medium <- terra::sel(tsi_list[[4]][[c(430:430)]])
# small  <- terra::sel(tsi_list[[4]][[c(430:430)]])

# used to make raster stack (c() stack does not work due to extents not matching)
#roi_list <- c(watt_GL, large, medium,small)


#save roi sample points raster stack

# writeRaster(watt_GL, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/watt_GL.tif",
#           filetype='GTiff', overwrite=TRUE)
# writeRaster(large, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/large.tif",
#           filetype='GTiff', overwrite=TRUE)
# writeRaster(medium, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/medium.tif",
#           filetype='GTiff', overwrite=TRUE)
# writeRaster(small, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/small.tif",
#           filetype='GTiff', overwrite=TRUE)

#load sampled points
watt_GL <- rast("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/watt_GL.tif")
large <- rast("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/large.tif")
medium <- rast("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/medium.tif")
small <- rast("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/roi_sample_locations/small.tif")

# used to make raster stack (c() stack does not work due to extents not matching)
#roi_list <- c(watt_GL, large, medium,small)
roi_list <- list(watt_GL, large, medium,small)
 
# Convert each ROI to df
# 
# watt_GL_df <- watt_GL  |>  as.data.frame(xy=T) |> select(x,y)    |> mutate(type="watt_GL")
# large_df   <- large    |>  as.data.frame(xy=T) |> select(x,y)    |> mutate(type="large")
# medium_df  <- medium   |>  as.data.frame(xy=T) |> select(x,y)    |> mutate(type="medium")
# small_df   <- small    |>  as.data.frame(xy=T) |> select(x,y)    |> mutate(type="small")

# Combine locations into single df
# sampled_points <- rbind(watt_GL_df, large_df, medium_df, small_df)

#save sampled points
#saveRDS(sampled_points, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/sampled_points.rds")

# load sampled points
sampled_points <- readRDS("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/sampled_points.rds")

Plots of ROIs

sampled location (little red highlighed areas)

knitr::include_graphics("~/Documents/R/HU_SHK/tsi_1984-2021/ROI_location_map.png")

Detailed plots of each ROI

lapply(seq_along(roi_list),
       function(i)
         plot(tsi_list[[4]][[c(430:430)]],
              legend=T, 
              col = viridis(n=100,option="D"),
              range=c(0,10000),
              ext=ext(roi_list[[i]])*10,
             plg=list(title= paste(roi_names[i], "\n", names(tsi_list[[i]][[600:600]]))
                                 ))
)

compairison of TSI cases for each location

lapply(seq_along(roi_list),
       function(i)
         plot(tsi_stack,
              legend=T, 
              col = viridis(n=100,option="D"),
              ext=(ext(roi_list[[i]])*10),
               plg=list(title= paste(names(tsi_list[i]), "\n", roi_names[i])),
              range=c(0,10000))
       
)

For the displayed date (July 2013), all plots seem to show quite congruous results. As expected the un-interpolated time series shows more data gaps, but in general the visual correspondence of reflectance values look good.

TSI visualization

data manipulation and subsample creation

Due to C memory usage limit issues, the conversion from spatraster to spatial data frame had to be split into small individual processing chunks, that are then subsequently joined to form the single df. First unique start and stop limits for the processing loop are defined (as cases have heterogenious lengths depending on number of NA layers and date window size). Following this the range of the processing checks is set in a list. Then the spat raster object of each case are converted into spatial df´s, via the defined subset processing chunks. Finally, all the individual chunks of a case are joined into a single spatial df object

… works but is not a pretty workaround…

tsi_df   <- list("base_case"=c(), "day_16_sigma_8_16_32"=c(), "day_16_sigma_4_8_16"=c(), "day_16_sigma_8_8_8_16_16_16_64"=c(), "day_8_sigma_8_16_32"=c())

start <- list()
stop <- list()
for (i in 1:length(names(tsi_list))) {
  # create processing chunks 
  start[[i]] <- seq(from = 1, to = length(names(tsi_list[[i]]))-20, by = 20)
  stop[[i]] <- seq(from = 20, to = length(names(tsi_list[[i]])), by = 20)
  
}



# create a unique length of processing chuncks based on the number of observations in each scenario
processing_chunks <- start[[1]][c(1:length(start))] # initial list
# loop though each case
for (j in 1:length(start)) {
  # loop through each start and stop combination
  for (i in 1:length(start[[j]])) {
    
    processing_chunks[[j]][i] <- list(c(start[[j]][i], stop[[j]][i]) )
    
  }
  
}

# Convert spat raster into df, via subset processing chunks

### uncomment to execute
# lapply(seq_along(tsi_list),
#        function(i)
#
#          for (j in c(1:length(processing_chunks[[i]]))) {
#
#            tsi_df[[i]][[j]] <<- terra::as.data.frame(tsi_list[[i]][[c(processing_chunks[[i]][[j]][1]:
#                                                                         processing_chunks[[i]][[j]][2])]],
#                                                      xy = TRUE, na.rm=F) |>
#              # filter out rows that only contain na
#              filter_at(vars(-x,-y), any_vars(! is.na(.))) |>
#
#              left_join(sampled_points) |>
#               filter(!is.na(type))
#          }
# )


# create join for list into df
# tsi_df_joined <- list()
#  lapply(seq_along(tsi_df),
#         function(i)
#          tsi_df_joined[[i]] <<-  reduce(tsi_df[[i]], full_join) |>
#           group_by(type) |>
#           sample_n(1)
#  )

#save tsi df
#saveRDS(tsi_df_joined, file = "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/tsi_df_joined.rds")

# load tsi df
tsi_df_joined <- readRDS("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/R_objects/tsi_df_joined.rds")

# create long version of df with unique and formatted location and date cols
tsi_df_long <- list()
lapply(seq_along(tsi_list),
       function(i) 
         tsi_df_long[[i]] <<- tsi_df_joined[[i]] |> 
         pivot_longer(cols = c(-x,-y,-type), names_to = "timestamp") |> 
         mutate(xy= paste0(x,"_",y)) |> 
         mutate(date = as.Date(timestamp, format="%Y%m%d"),
                year = lubridate::year(date),
                value = if_else(value<0,0,as.numeric(value)),
                case  = names(tsi_list[i]))
       
)

plot faceted time series 1984:2000

lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_df_long[[i]] %>%
         filter(year<2001) %>%
         
         ggplot( aes(date, value, color=type, group=xy)) +
         geom_line(aes(group=type)) +
         geom_point(size=0.5) +
         theme_classic() +
         scale_colour_viridis_d(option = "D") +
         scale_x_date(date_breaks = "1 month", date_labels = "%m")+
         theme_minimal() +
         theme(axis.text.x = element_text(angle = 45,  hjust=1),
               legend.position="bottom"
         ) +
         ggtitle(cases[i]) +
         facet_wrap(~year, scales = "free_x")
       
)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

plot faceted time series 2001:2021

lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_df_long[[i]] %>%
         filter(year>2001) %>%
         
         ggplot( aes(date, value, color=type, group=xy)) +
         geom_line(aes(group=type)) +
         geom_point(size=0.5) +
         theme_classic() +
         scale_colour_viridis_d(option = "D") +
         theme_minimal() +
         scale_x_date(date_breaks = "1 month", date_labels = "%m")+
         theme(axis.text.x = element_text(angle = 45,  hjust=1),
               legend.position="bottom"
         ) +
         ggtitle(cases[i]) +
         facet_wrap(~year, scales = "free_x")
       
)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

lapply(seq_along(tsi_list),
       function(i) 
         ggplot(tsi_df_long[[i]], aes(date, value)) +
         geom_smooth() +
         theme_minimal() +
         scale_x_date(date_breaks = "1 year", date_labels = "%d-%m-%y")+
         theme(axis.text.x = element_text(angle = 45 ,  hjust=1)) 
)
lapply(seq_along(tsi_list),
       function(i) 
         
         tsi_df_long[[i]] %>%
         
         
         ggplot( aes(date, value, group = year)) +
         geom_boxplot() +
         theme_minimal() +
         scale_x_date( date_labels = "%d-%m-%y")+
         theme(axis.text.x = element_text(angle = 45, hjust=1)) 
)

direct location and case compairison

tsi_df_long_combined <- reduce(tsi_df_long, full_join)
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
## Joining, by = c("x", "y", "type", "timestamp", "value", "xy", "date", "year",
## "case")
ggplot(tsi_df_long_combined, aes(date, value, color=case, group = xy)) +
  geom_line(aes(group=case), alpha=0.5) +
  geom_point(size=0.5, alpha=0.9) +
  geom_smooth(color="black") +
  theme_minimal() +
  scale_colour_viridis_d(option = "D") +
  scale_x_date(date_breaks = "4 year", date_labels = "%Y")+
  theme(axis.text.x = element_text(angle = 45,  hjust=1),
        legend.position="bottom") +
  #  ggtitle(cases[i]) +
  facet_grid(rows = vars(type), cols = vars(case), scales="free_y") 
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#facet_grid( ~year)
#facet_grid(vars(drv), vars(cyl))

compairson of Base case values and interpolated TSI outputs

year_chunks <- round(seq(1984, 2021, length.out = 5))
lapply(seq_along(c(1:4)),
       function(i) 

tsi_df_long_combined %>%
  filter(year > year_chunks[i],
         year < year_chunks[i+1],
         case != "base_case") %>%
  
  ggplot(aes(date, value, color=case, group = xy)) +
  geom_line(aes(group=case)) +
  geom_point(size=0.5) +
  geom_point(data=tsi_df_long_combined |> filter(case =="base_case",
                                                 year > year_chunks[i],
                                                 year < year_chunks[i+1]),
             aes(date, value), 
             color ="red", alpha=0.5) +
  theme_minimal() +
  scale_colour_viridis_d(option = "D") +
  scale_x_date(date_breaks = "1 month", date_labels = "%m")+
  theme(axis.text.x = element_text(angle = 45,  hjust=1),
        legend.position="bottom") +
  facet_grid(rows = vars(type), cols = vars(year), scales="free_x") 

#  facet_wrap(~year, scales = "free_x") 

)
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

plot time series with formatted dates

ggplot(tsi_df_long[[i]], aes(date, value, group=xy, color=xy)) +
  geom_line() +
  theme_classic() +
  scale_colour_viridis_d(option = "D") +
  scale_x_date(date_breaks = "1 month", date_labels = "%m-%y")+
  theme(axis.text.x = element_text(angle = 45,  hjust=1),
        legend.position="none")

ggplot(tsi_df_long[[i]], aes(date, value)) +
  geom_smooth() +
  theme_minimal() +
  scale_x_date(date_breaks = "1 year", date_labels = "%d-%m-%y")+
  theme(axis.text.x = element_text(angle = 45 ,  hjust=1)) 

ggplot(tsi_df_long[[i]], aes(date, value, group = timestamp)) +
  geom_boxplot() +
  theme_minimal() +
  scale_x_date( date_labels = "%d-%m-%y")+
  theme(axis.text.x = element_text(angle = 45, hjust=1)) 

ggplot + tidyterra vis option (nice but runs slow)

ggplot() +
  geom_spatraster(data = tsi[[c(4:8)]]) +
  facet_wrap(~lyr, ncol = 2) +
  scale_fill_whitebox_c(
    palette = "viridi",
    # labels = scales::label_number(suffix = "º")
  ) +
  labs(fill = "NDVI") +
  theme_classic()
## gif image creation

for (i in c(1:100
            #length(names(tsi))
)) {
  
  # Step 1: Call the pdf command to start the plot
  png(file = paste0("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/day_16_sigma_8_16_32/plot_", names(tsi_list[[1]])[[i]], ".png"), width = 265, height = 125, units='mm', res = 100) #,   # The directory you want to save the file in
  # width = 12, # The width of the plot in inches
  #height = 12) # The height of the plot in inches
  plot(tsi_list[[2]], legend=T, col = viridis(n=100,option="D"),
       range=c(0,10000),
       plg=list( title = paste(as.Date(names(tsi_list[[1]][[i]]), format="%Y%m%d")), title.cex=1.25)
  )
  
  
  # Step 3: Run dev.off() to create the file
  dev.off()
}

TSI gif

knitr::include_graphics("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/day_16_sigma_8_16_32/ts_16_8_16_32.gif")

# old attempted at looped df processing

processing_chunks <- list(c(1,10),c(11,20),c(21,30),c(31,40),c(41,50))


processing_chunks <- c(1,100,200,300,400,500,600,700,800)
processing_chunks <- c(1,10,20,30) # short version for testing

tsi_df   <- list("base_case"=c(), "day_16_sigma_8_16_32"=c(), "day_16_sigma_4_8_16"=c(), "day_8_sigma_8_16_32"=c())

for (j in c(1:3
            #(length(processing_chunks)-1)
)) {
  
  lapply(seq_along(tsi_list),
         function(i) 
           
           tsi_df[[i]][[j]] <<- terra::as.data.frame(tsi_list[[i]][[c(processing_chunks[[j]][1]:
                                                                        processing_chunks[[j]][2])]], 
                                                     xy = TRUE, na.rm=F) |>
           # filter out rows that only contain na
           filter_at(vars(-x,-y), any_vars(! is.na(.))) |> 
           # random subset for data vis
           slice_sample(n=15)
  )
  
}

foo <- list()
lapply(seq_along(tsi_df),
       function(i) 
         
         foo[[i]] <<-  reduce(tsi_df[[i]], left_join)
       
)
### test for R tile merging 

link_a <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0057_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

link_b <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/day_16_sigma_8_16_32/X0054_Y0052/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"

a <- rast(link_a)
b <- rast(link_b)

ab <- merge(a,b)

plot(ab[[c(10:13)]], legend=F, col = viridis(n=100,option="D"),
     range=c(0,10000)
)

object.size((ab))

writeRaster(ab, "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/test_mosaic_r.tif",
            # overwrite=TRUE,
            datatype = "INT4S",
            filetype='GTiff')


link_ab <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/test_mosaic_r.tif"

ab <- rast(link_ab)
plot(ab[[c(10:13)]], legend=F, col = viridis(n=100,option="D"),
     range=c(0,10000)
)